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A new approach to turbulence modeling 


By B. Perot 1 AND P, Moin 2 


A new approach to Reynolds averaged turbulence modeling is proposed which has 
a computational cost comparable to two equation models but a predictive capability 
approaching that of Reynolds stress transport models. This approach isolates the 
crucial information contained within the Reynolds stress tensor, and solves trans- 
port equations only for a set of “reduced” variables. In this work, direct numerical 
simulation (DNS) data is used to analyze the nature of these newly proposed tur- 
bulence quantities and the source terms which appear in their respective transport 
equations. The physical relevance of these quantities is discussed and some initial 
modeling results for turbulent channel flow are presented. 


1. Introduction 

1.1 Background 

Two equation turbulence models, such as the k/e model and its variants, are 
widely used for industrial computations of complex flows. The inadequacies of these 
models are well known, but they continue to retain favor because they are robust 
and inexpensive to implement. The primary weakness of standard two equation 
models is the Boussinesq eddy viscosity hypothesis: this constitutive relationship is 
often questionable in complex flows. Algebraic Reynolds stress models (or non-linear 
eddy viscosity models) assume a more complex (nonlinear) constitutive relation for 
the Reynolds stresses. These models are derived from the equilibrium form of the 
full Reynolds stress transport equations. While they can significantly improve the 
model performance under some conditions, they also tend to be less robust and 
usually require more iterations to converge (Speziale, 1994). The work of Lund &: 
Novikov (1992) on LES subgrid closure suggests that even in their most general 
form, non-linear eddy viscosity models are fundamentally incapable of completely 
representing the Reynolds stresses. Industrial interest in using full second moment 
closures (the Reynolds stress transport equations) is hampered by the fact that 
these equations are much more expensive to compute, converge slowly, and are 
susceptible to numerical instability. 

In this work, a turbulence model is explored which does not require an assumed 
constitutive relation for the Reynolds stresses and may be considerably cheaper to 
compute than standard second moment closures. This approach is made possible 
by abandoning the Reynolds stresses as the primary turbulence quantity of interest. 
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The averaged Navier-Stokes equations only require the divergence of the Reynolds 
stress tensor, hence the Reynolds stress tensor carries twice as much information as 
required by the mean flow. Moving to a minimal set of turbulence variables reduces 
the overall work by roughly half, but introduces a set of new turbulence variables, 
which at this time are poorly understood. This project attempts to use DNS data 
to better understand these new turbulence variables and their exact and modeled 
transport equations. 

1.2 Formulation 

The averaged Navier-Stokes equations take the following form for incompressible, 
constant-property, isothermal flow: 


V ■ u = 0 (la) 

du 

— + u • Vu = -Vp + uV • S - V * R (lft) 

where u is the mean velocity, p is the mean pressure, v is the kinematic viscosity, 
S = Vu + (Vu) T is twice the rate-of-strain tensor, and R is the Reynolds stress 
tensor. The evolution of the Reynolds stress tensor is given by: 

■^ + uVR= r/V 2 R + P- e + n- V- T-[Vq + ( Vq) r ] (2) 


where P is the production term, e is the (homogeneous) dissipation rate tensor, 
n is the pressure-strain tensor, T is the velocity triple-correlation, and q is the 
velocity-pressure correlation. The last four source terms on the right-hand side 
must be modeled in order to close the system. The production term P is exactly 
represented in terms of the Reynolds stresses and the mean velocity gradients. This 
is the standard description of the source terms, but it is by no means unique and 
there are numerous other arrangements. 

Note that turbulence effects in the mean momentum equation can be represented 
by a body force f = V R. One could construct transport equations for this body 
force (which has been suggested by Wu et a/., 1996), but mean momentum would no 
longer be simply conserved. To guarantee momentum conservation, the body force 
is decomposed using Helmholtz decomposition, into its solenodal and dilatational 
parts, f = V0 + V x ip. A constraint (or gauge) must be imposed on ip to make the 
decomposition unique. In this work we take V • ip = 0. With this choice of gauge, 
the relationship between <f> and ip and the Reynolds stress tensor is given by, 

V 2 </> = V • (V • R) (3a) 

V 2 v> = -V X (V • R) (36) 

Note that the choice of gauge influences the value of ip , but does not affect how ip 
influences the mean flow. 
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Using these relationships, transport equations for <j> and ip can be derived from 
the Reynolds stress transport equations. 


d 4- + U • V<f> = i /VV - 2V • q - V' 2 V •V-[e-n + V- T-P] + V~ 2 S< 
at 


^ + u • VV> = i/V 2 V> + V x q + V" 2 V xV-[e — II + V- T — P] + V 2 S^ 
at 


(4a) 

(46) 


These equations contain extra production-like source terms S# and which contain 
mean velocity gradients. Note that the production term is not an explicit function 
of (j> and tj) (except under limited circumstances) and, in general, must be modeled. 
The inverse Laplacian V“ 2 that appears in these equations can be thought of as an 
integral operator. 


2. Theoretical analysis 

2.1 Turbulent pressure 

Taking the divergence of Eq. (lb) (the mean momentum equation) gives the 
classic Poisson equation for pressure, 

V 2 p = —V * (u * Vu) — V • (V • R) (5) 

Since this is a linear equation, the pressure can be split conceptually into two terms: 
one can think of the mean pressure as being a sum of a mean flow pressure due to 
the first term on the right-hand side, 

V 2 P mefln = -V-(u-Vu) (6a) 

and a turbulent pressure due to the second term on the right-hand side, 

V 2 Pturb = -V - (V • R) (66) 

Given the definition of <f) and assuming that <f> is zero when there is no turbulence, 
then it is clear that <f> = — P tU rb • For this reason, <j> will be referred to as the 
turbulent pressure. This quantity is added to the mean pressure in the averaged 
momentum equation, which results in P me an = p + <f> being the effective pressure 
for the averaged equations. The quantity P mta n tends to vary more smoothly than 
p, which aids the numerical solution of these equations. 

For turbulent flows with a single inhomogeneous direction, the turbulent pressure 
can be directly related to the Reynolds stresses. In this limit Eq. (3a) becomes 
<j> 22 = #22,22 where X2 is the direction of inhomogeneity. This indicates that </> = 
#22 for these types of flows. Note that #22 is positive semi-definite, so <f> is always 
greater than or equal to zero in this situation. Positive <f> is consistent with the 
picture of turbulence as a collection of random vortices (with lower pressure cores) 
embedded in the mean flow. It is not clear what the conditions for a negative 
turbulent pressure would be, if this condition is indeed possible. 
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2.2 Turbulent vorticity 

To understand the role of ip it is instructive to look again at turbulent flows that 
have a single inhomogeneous direction. Under this restriction Eq. (3b) becomes 
ipi t 22 = — £ i 2 ,22 where X 2 is the direction of inhomogeneity. If ip goes to zero 
when there is no turbulence then ipi — — ti 2 kRk 2 - (or i/q — — i? 32 , ip 2 = 0 and 
if ? 3 = P\2 )* These are the off diagonal, or shear stress components of the Reynolds 
stress tensor. 

For two-dimensional mean flows with two inhomogeneous flow directions, only 
the third component of ip is non-zero, and Eq. (3b) becomes 

V ? 3,ll + 03,22 = # 12,22 ” -^ 12, 11 + {R \ 1 ~ # 22), 12 (7) 

Since ip is responsible for vorticity generation, it is appropriate that it be aligned 
with the vorticity in two-dimensional flows. As a first level of approximation, it is 
not unreasonable to think of ip as representing the average vorticity of a collection 
of random vortices making up the turbulence, and therefore ip will be referred to 
as the turbulent vorticity. 

For two-dimensional flows with a single inhomogeneous direction xp$ = R\i- 
Note how the components of ip reflect the dimensionality of the problem, while the 
mathematical expressions for these components reflects the degree of inhomogeneity. 

2.3 Relationship with the eddy viscosity hypothesis 
The linear eddy viscosity hypothesis for incompressible flows takes the form, 

R= -i/ T (Vu + (Vu) T ) + hi (8) 

o 

where ut is the eddy viscosity, I is the identity matrix, and k is one half the trace 
of the Reynolds stress tensor. 

Taking the divergence of Eq. ( 8 ) and rearranging terms gives, 

2 

f = V • R = V( — 2u • Vut) + Vx (i/^V x u) + 2u • V( Wr ). (9) 

o 

If the eddy viscosity varies relatively slowly, as is usually the case, then the very last 
term (involving the second derivative of the eddy viscosity) will be small and can be 
neglected. Under these circumstances the linear eddy viscosity model is equivalent 
to the following model, 

2 

<P = -k — 2u • Vut (10a) 

o 

ip = i /yV x u. (106) 

So to a first approximation the turbulent vorticity, ip should be roughly equal to the 
mean vorticity, times a positive eddy viscosity; and the turbulent pressure should 
be roughly equal to two thirds of the turbulent kinetic energy. These results are 
entirely consistent with the findings of the previous subsections. 
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FIGURE 1. Contours of turbulent pressure (0) and negative turbulent vorticity 
(— 0) for the separating boundary layer of Na &: Moin. 

3. Computational results 

Equations (3a) and (3b), relating the turbulent pressure* and turbulent, vorticity 
to the Reynolds stresses, were used to calculate 0 and 0 from DNS data for two 
relatively complex two-dimensional turbulent flows: a separating boundary layer 
(Na & Moin, 1996) and flow over a backward facing step (Le & Moin, 1995). The 
purpose was to assess the behavior of these turbulence quantities in practical tur- 
bulent situations, and to provide a database of these quantities for later comparison 
with turbulence models. 


3.1 Separated boundary layer 

The values of 0 and —03 are shown in Fig. 1. As mentioned previously, for two- 
dimensional flows only the third component of 0 is nonzero. The flow moves from 
left to right, separates just before the midpoint of the computational domain, and 
then reattaches before the exit. The contours are the same for both quantities and 
range from -0.0004/7^, to O.OIL^, where Uoo is the inlet free- stream velocity. 

Both the turbulent pressure and turbulent vorticity magnitudes increase in the 
separating shear layer and the reattachment, zone. In addition, both quantities 
become slightly negative in the region just in front (to the left) of the separating 
shear layer, and show some “elliptic” (long range decay) effects at the top of the 
separation bubble. There is some speculation at this time that these effects could 
be numerical, but there is also some reason to believe that they are a legitimate 
result of the elliptic operators which define these variables. Changes in the far- 
field boundary condition (from zero value to zero normal gradient) had no visibly 
perceptible effect on the values of 0 and 0 3 . 

The visual observation that 0 and — 0 3 are roughly proportional is analogous to 
the observation that 0.3A* « ~Rn (originally developed by Townsend, 1956, and 
successfully used in the turbulence model of Bradshaw, Ferriss &, Atwell, 1967). 
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FIGURE 2. Contours of the normal Reynolds stress (R 2 2 ) and negative turbulent 
shear stress ( — R\ 2 ) for the separating boundary layer of Na & Moin. 


It is also consistent with the (first order) notion of turbulence as a collection of 
embedded vortices, with —<t> representing the average vortex core pressure and ip 
representing the average vortex strength. 

In the case of a single inhomogeneous direction, <p ~ R 22 and 0 3 — R\ 2 . It is 
instructive therefore to compare the results shown in Fig. 1 with the R 22 and —R\ 2 
components of the Reynolds stress tensor, shown in Fig. 2. The magnitudes of 
the contours in Fig. 2. are the same as Fig. 1. This comparison clearly shows the 
additional effects that result from inhomogeneity in the streamwise direction. The 
leading and trailing boundary layers (which have very little streamwise inhomogene- 
ity) are almost identical. However, the magnitudes of the turbulent pressure and 
turbulent vorticity are enhanced in the separated shear layer due to the streamwise 
inhomogeneity. 

3.2 Backward facing step 

Computations of <j> and — for the backward facing step are shown in Fig. 3. The 
flow is from left to right, and there is an initial (unphysical) transient at the inflow 
as the inflow boundary condition becomes Navier-Stokes turbulence. The boundary 
layer leading up to the backstep has moderate levels of the turbulent pressure and 
turbulent vorticity (which closely agree with the values of R 22 and — i? i2 in that 
region). As with the separating boundary layer, the turbulent pressure and turbu- 
lent vorticity increase significantly in the separated shear layer and reattachment 
zone. There is an area of slight positive turbulent pressure and negative turbulent 
vorticity in the far field (about one step height) above the backstep corner. This 
may or may not be a numerical artifact, and is discussed in the next section. 

3.3 Ellipticity 

Identifying the exact nature of the ellipticity of these new turbulence quantities is 
important to understanding their overall behavior and how they should be modeled. 


A new approach to turbulence modeling 


41 



PHI 



PSI 

Figure 3. Contours of turbulent pressure (< p ) and negative turbulent vorticity 
(—ip) for the backward facing step of Le &: Moin. 

When rewritten, Eqs. (3a) and (3b) become, 

<t> = V" 2 V.(V*R) (1 la) 

tp = -V- 2 V x (V • R) (116) 

These are elliptic, but order one, operators on the Reynolds stress tensor. As demon- 
strated in §2, when there is only a single inhomogeneous direction, these operators 
simply lead to various Reynolds stress components. Under these conditions they do 
not produce “action at a distance” or long range effects normally associated with 
elliptic (Poisson or Helmholtz) operators. 

For two and three inhomogeneous directions, it is still not clear whether these 
operators produce long range effects. There are certainly some situations in which 
they do not. One example is when the Reynolds stress tensor can be represented 
in the following form (somewhat reminiscent of the linear eddy viscosity relation) 
Rtj = s6 t j + Vij + vj yl , where s is some scalar and v is a vector. If this is the case 
then, <f> = s + 2V ■ v and ip — — V x v, and there are no long range (“elliptic”) 
effects. 

In fact, the presence of long range effects in <p and ip is somewhat unsettling. It 
would suggest that these turbulence quantities can exist in regions where there is 
no Reynolds stress. Since V * R = + V x ip, this would imply that a precise 

cancellation of these long range effects must occur in regions where the Reynolds 
stresses are small or negligible. While the results presented in Fig. 1 and Fig. 3 
seem to show that long range elliptic effects do indeed take place, they could also 
be a numerical artifact. The numerical solution of Eqs. (3a) and (3b) requires 
double differentiation of the DNS data; this produces compact Poisson equation 
source terms that are only marginally resolved by the mesh. It is our current 
conjecture that these operators are actually local in nature and only serve to “mix” 
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FIGURE 4. Budget of the <f> transport equation at a station roughly half way 
through the recirculation bubble of the backward facing step ( xj h - 4.0). - dis- 
sipation or diffusion; velocity pressure-gradient; triple correlation term; 

production (positive) or convection. 

various components of the Reynolds stress tensor. It is also conjectured from these 
computational results that the turbulent pressure is a positive semi-definite quantity. 

Note that the ellipticity discussed here is not the same as an ellipticity in the 
governing evolution equations for these quantities. An elliptic term in the evolution 
equations is both physical and desirable (see Durbin, 1993). Such a term mimics 
long range pressure effects known to occur in the exact source terms. The exact 
evolution equations for <t> and Vs described below, have just this elliptic property. 

3.4 Turbulent pressure evolution 

Considerable insight can be obtained about the evolution of the turbulent pres- 
sure by considering the case of a single inhomogeneous direction. It has been shown 
that under these circumstances <t> = i?22i so the evolution is identical with the 
Reynolds stress transport equation for the normal Reynolds stress, i? 22 . For the 
case of turbulent channel flow (Mansour et al , 1988), the R 2 2 evolution is domi- 
nated by a balance between dissipation and pressure-strain, with somewhat smaller 
contributions from turbulent transport and viscous diffusion. There is considerable 
interest in determining if these same trends continue for <j> evolution in more com- 
plex situations, since the ultimate goal is to construct a modeled evolution equation 
for this quantity. 

Figure 4 shows the terms in the exact <j> evolution equation for flow over a back- 
ward facing step, at a station roughly in the middle of the recirculation bubble. 
These terms were calculated in the same manner as the turbulent pressure. Both 
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FIGURE 5. Budget of the 03 transport equation at a station roughly half way 
through the recirculation bubble of the backward facing step: see Fig. 4 for caption. 


the detached shear layer and the backward moving boundary layer are visible in the 
statistics. In the shear layer, the expected dominance of dissipation and pressure- 
terms (presumably dominated by pressure-strain) is evident. In the recirculating 
boundary layer, turbulent transport and pressure-terms (probably dominated by 
pressure transport) are dominant. It is interesting to note that the production term 
dominates in the middle of the recirculation bubble. The fact that some of these 
source terms are not exactly zero at roughly two step heights away from the bottom 
wall is thought to be a numerical artifact similar to those found when calculating 
(f> and 0. Some of the curves have an erratic nature due to the lack of statisti- 
cal samples. This phenomena is also present in the (unsmoothed) Reynolds stress 
transport equation budgets presented in Le & Moin, 1993. 

3.5 Turbulent vorticity evolution 

As with the turbulent pressure, it is useful to consider the case of a single inhomo- 
geneous direction when analyzing the evolution of the turbulent vorticity. Under 
these circumstances 03 evolves identically to the Reynolds shear stress, i?i2 * In 
turbulent channel flow, the R 12 evolution is dominated by a balance between pro- 
duction and pressure-strain, with somewhat smaller contributions from turbulent 
and pressure transport. This trend continues in the 03 evolution equation, which 
is shown in Fig. 5., for the backward facing step at a cross section roughly halfway 
through the recirculization bubble ( x/h = 4.0). The small value of the dissipation 
is consistent with the fact that isotropic source terms can be shown not contribute 
to the evolution of 0. 
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4. Modeling 

4 A Formulation 

An initial proposal for modeled transport equations for the turbulent pressure 
and turbulent vorticity are, 


d<p 

dt 


+ u- V0 = V-(u + ur)¥<p 


dip 


-5— + u ■ Vip = V • (u + uj)V ip 
at 


-(NGM?)-*© 


t P • Ip 
15 u + ur 


(12 a) 
(126) 


where, C ^ = 0.09, y is the normal distance to the wall, the time-scale is given by T = 
(u + u T )/<p, and the eddy viscosity is given by Ur — |V>|/M- Dissipation (and some 
redistribution) is modeled as an exponential decay process (roughly corresponding to 
Rotta’s, low Reynolds number dissipation model). Turbulent and pressure transport 
are collectively modeled as enhanced diffusive transport. Production and energy 
redistribution are proportional to the turbulence pressure times the mean vorticity 
for the turbulent vorticity, and are proportional to the square of the turbulent 
vorticity magnitude for the turbulent pressure. High Reynolds number constants are 
determined so that <p = ~k at high Reynolds numbers. The low Reynolds number 
constants (which appear with a u) are set to obtain exact asymptotic behavior and 
good agreement with the channel flow simulations of the next section. 

Note that both <f> and ip have the same units. An extra turbulent scale is currently 
defined by using the mean flow timescale \u\ to define the eddy viscosity. The 
solution of an additional scale transport equation (such as e) would remedy a number 
of potential problems with the current model. It could eliminate the singularity in 
the eddy viscosity at zero vorticity, remove any explicit references to the wall normal 
distance, and allow better decay rates for homogeneous isotropic turbulence. The 
disadvantage of this approach (which will be tested in the future) is the added 
computational cost and additional empiricism. 


4>2 Channel flow simulations 

The model equations (12a and 12b) were solved in conjunction with mean flow 
equations for fully developed channel flow at Re r of 180 and 395. Since there is only 
one inhomogeneous direction, the turbulent pressure is proportional to the normal 
Reynolds stress, and is proportional to the turbulent shear stress. Comparisons 
of the model predictions and the DNS data of Kim, Moin, & Moser (1987), are 
shown in Fig. 6. 

When a turbulent channel flow is suddenly perturbed by a spanwise pressure gra- 
dient, the flow suddenly becomes three dimensional and the turbulence intensities 
first drop before increasing due to the increased total shear (Moin et al , 1990). 
Durbin (1993) modeled this effect by adding a term to the dissipation equation 
which increases the dissipation in these three-dimensional flows. The same quali- 
tative effect can be obtained by defining the eddy viscosity in the proposed model 

as ur = I n two-dimensional flows this is identical to the previous definition. 
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FIGURE 6. Model results (solid lines) and DNS data (circles) for turbulent channel 
flow. (Re T = 180) 

However, in three-dimensional flows, the orientation of ip will lag u>, and the eddy 
viscosity will drop initially. A smaller eddy viscosity leads to a smaller timescale 
and increased dissipation. Unfortunately, the magnitude of this effect is severely 
underestimated in the present model, and a scale equation (and a correction like 
Durbin’s) may be required to model this effect accurately. 

5. Conclusions 

This work proposes abandoning the Reynolds stresses as primary turbulence 
quantities in favor of a reduced set of turbulence variables, namely the turbulent 
pressure cp, and the turbulent vorticity ip. The advantage of moving to these al- 
ternative variables is the ability to simulate turbulent flows with the accuracy of a 
Reynolds stress transport model (i.e. with no assumed constitutive relations), but 
at a significantly reduced cost and simplified model complexity. As the names im- 
ply, these quantities are not simply mathematical constructs formulated to replace 
the Reynolds stress tensor. They are physically relevant quantities. 

At first glance the operators which relate <f> and ip to the Reynolds stress tensor 
suggest the possibility of ellipticity or action at a distance. However, we have 
shown that under a number of different circumstances this does not happen, and 
conjecture that it may never happen. The physical relevance of these quantities 
would be complicated if they were finite when there was no turbulence (Reynolds 
stresses). A proof to this effect may also prove our second conjecture, that <p is a 
positive definite quantity. 

The budgets for the transport equations of these new variables indicated that 
the extra production terms were not significant, and that these equations could be 
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modeled analogously to the Reynolds stress transport equations. An initial model 
was constructed for these equations using basic modeling constructs which showed 
good results for turbulent channel flow. It is likely, that for this shearing flow, the 
turbulent timescale is well represented by the mean flow vorticity. However, for 
more complex situations, it is likely that an additional scale equation (such as an e 
equation) will be required. 
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